#include "libmri.h"

void calib3d_cheev (complex *&V, int& nv,
                    complex *calib, int *csize, int *ksize, float tau)
{
	int m, n;

	calib3d_matrix_size (m, n, csize, ksize);

	complex *A = new complex [m*n];
	complex *AHA = new complex [n*n];
	float *L = new float [n];
	int info;

	calib3d_matrix ((cplxf*) A, (cplxf*) calib, csize, ksize);

	complex  one = {1.0f, 0.0f},
	        zero = {0.0f, 0.0f};

	for (int i = 0; i < n*n; ++i)
		AHA[i].real = AHA[i].imag = 0.0f;

	cgemm ('C', 'N', n, n, m, &one, A, m, A, m, &zero, AHA, n);
	cheev ('V', 'U', n, AHA, n, L, &info);

	float maxL = fabsf (L[0]);
	
	for (int i = 1; i < n ; ++i)
		if (fabsf (L[i]) > maxL)
			maxL = fabsf (L[i]);

	nv = 0; int *vi = new int[n];
	for (int i = 0; i < n ; ++i)
		if (fabsf (L[i]) >= tau*maxL)
			 vi[nv++] = i;

	V = new complex [nv*n];

	for (int j = 0; j < nv; ++j)
		for (int i = 0; i < n; ++i)
			V [i + n*j] = AHA [i + n*vi[j]];

	delete [] vi;
	delete [] A;
	delete [] AHA;
	delete [] L;
}

